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ABSTRACT 


Approximate  algorithms  have  long  been  the  only  available  methods 
for  generating  Poisson  random  variates  when  the  mean  is  large.  A  new 
algorithm  is  developed  which  is  exact,  has  execution  time  which  is 
insensitive  to  the  value  of  the  mean,  and  is  valid  whenever  the  mean  is 
greater  than  ten.  This  algorithm  is  compared  to  the  three  other 
algorithms  which  have  been  developed  recently  for  generating  Poisson 
variates  when  the  mean  is  large.  Criteria  used  are  set-up  time, 
marginal  execution  time,  memory  requirements,  and  lines  of  code.  New 
simple  tight  bounds  on  Poisson  probabilities  contribute  to  the  speed  of 
the  algorithm,  but  are  useful  in  a  more  general  context.  In  addition,  a 
survey  of  Poisson  variate  generation  algorithms  is  given. 


1.  INTRODUCTION 


We  consider  algorithms  for  generating  random  variates  from  the 
Poisson  mass  function 

fp(x)  =  e_li  »x/x!  x=0,1,2,... 

=  0  elsewhere, 

where  4  denotes  the  expected  value  of  the  random  variable  X.  In  Section 
2  existing  algorithms  for  Poisson  variate  generation  are  surveyed.  A 
new  algorithm,  PTPE,  is  developed  in  Section  3.  Computational  results 
are  shown  in  Section  4.  The  validity  of  PTPE  is  discussed  in  the 
Appendix . 

2.  LITERATURE  SURVEY 

Each  of  the  four  fundamental  approaches  to  variate  generation: 
inverse  transformation,  special  rties,  composition,  and 

acceptance/ rej ection,  (Schmeiser  C183)  has  b'.rn  used  as  the  basis  for 
existing  algorithms,  which  we  briefly  survey  here.  U(0,1>  is  used  to 
denote  the  uniform  distribution  over  the  unit  interval. 
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Probably  the  most  basic  approach  for  generating  random  variates  of 
any  kind  is  the  inverse  transformation. 

Algorithm  PINV 

1.  Generate  u  U(0,1),  set  x  *  0,  p  *  e  *i. 

2.  If  u  £  p ,  then  return  x. 

3.  Set  x  «-  x  +  1,  u  ♦  u-p,  p  *  py/x,  and  go  to  2. 

When  more  than  one  variate  is  to  be  generated  for  a  fixed  value  of  y, 

PINV  may  be  modified  to  save  the  initial  value  of  p  in  Step  1  and  the 

cumulative  probabilities  implicit  in  Step  3.  Either  way,  the  execution 

time  per  variate  increases  proportionally  with  y.  Fishman  Cl  1 □ 

1  /2 

developed  algorithm  PI F  which  executes  in  time  proportional  to  y  by 
performing  the  inverse  transformation  beginning  at  the  mode  and 
searching  either  increasingly  or  decreasingly  for  values  of  x.  To  begin 
the  search  at  the  mode,  both  the  cumulative  probability  p-CX  _<  y>  and 
probability  p-CX  =  y>  are  stored.  Fishman  stored  these  probabilities  for 
y=1,2, . . .,100  to  six  decimal  places,  but  the  size  and  accuracy  of  the 
table  could  easily  be  modified.  The  cumulative  probabilities  are 
calculated  recursively  as  in  PINV.  Snow  C20D  suggested  explicitly 
storing  the  cumulative  probabilities  and  using  binary  search  to 
determine  x.  Chen  and  Asau  C73  proposed  an  index  table  approach  (for 
a  y  discrete  distribution)  which  searches  the  cumulative  probabilities 
quickly  by  beginning  near  the  appropriate  value,  and  Atkinson  C5D 
included  an  algorithm  based  on  index  tables  in  his  computational 


results . 
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Special  properties  have  been  the  basis  for  several  Poisson 
algorithms.  The  best  known  and  simplest  is  based  on  the  exponential 
inter-event  times  of  the  homogeneous  Poisson  point  process. 

Algorithm  PMUL 

1 .  Set  x  <-0,  s  <-1,  p  «-  e  w. 

2.  Generate  u  v  11(0,1),  and  set  s  «-  su. 

3.  If  s  £  p,  then  return  x. 

4.  Set  x  ♦  x  +  1,  and  go  to  2. 

As  with  PINV  the  execution  time  increases  proportional ly  with  p  and 

storing  the  initial  value  of  p  for  future  use  is  reasonable  when  the 

value  of  y  does  not  change  each  time  a  variate  is  generated.  Note  that 

PINV  is  faster  than  PMUL  for  all  values  of  p  whenever  the  generation  of 

a  U (0,1 )  variate  requires  more  time  than  the  total  time  required  for  a 

division,  subtraction  and  a  storage  move.  The  authors  have  seen  no 

implementation  where  PMUL  was  faster  than  PINV. 

In  addition  to  the  inverse  transformation  methods,  composition  can 

be  used  as  the  basis  for  Poisson  algorithms.  Composition,  or 

probability  mixing,  is  used  in  variate  generation  by  returning  a  variate 

n 

from  f.(x)  with  probability  p^  when  f(x)  =  e  p^  f ^ < x) ,  where  n  may  be 

i=1 

finite  or  infinite  and  each  f. (x)  is  either  a  discrete  probability  mass 
function  or  a  density  function.  Let  I  be  a  Poisson  random  variable  with 
mean  x,  x  _>  p.  Then  a  binomial  random  variable,  arising  from  I  trials, 
each  having  probability  of  success  p/x,  has  a  Poisson  distribution  with 
mean  p.  The  proof  is  direct  by  noting 
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f  (x)  =  z  Ce" V/i!]C<  -J  )<y/X>x  (CX-y>/X)1_x  3  x=0,1,2,...  . 

The  advantage  to  this  composition  approach  is  that  e  u  does  not 

need  to  be  calculated  during  setup.  Usually  X  =  1  (Ahrens  and  Dieter 

[23  and  Fishman  Cl  ID)  with  the  resulting  algorithm  being  used  to  supply 
x  from  the  fractional  portion  of  y  when  y  is  not  integer.  A  reasonable 
implementation  for  y  £  1  "thins"  a  Poisson  variate  with  unit  mean. 
Using  PINV  to  generate  the  variate  with  a  mean  of  one  yields 

Algorithm  PTH  (y  <  1 ) 

1.  Generate  u  ^  U(0,1),  set  x  «-  0,  k  ♦  0,  p  «-  .367879441171. 

2.  If  u<  p,  then  return  x. 

3.  Set  k  k+1,  u  ♦  u-p,  p  «-  p/k.  Generate  v  U(0,1). 

If  v  <  y,  then  set  x  +•  x+1 .  Go  to  2. 

Fishman  Cl  1 3  gives  the  algorithm  in  a  form  assuming  the  cumulative 

probabilities  for  y  =  1  are  tabled.  A  similar  algorithm  can  be  created 

by  incrementing  x  in  Step  4  of  PWUL  with  probability  y  and  initializing 

-i 

p  ♦  .367879441171  =  e  in  Step  1.  The  idea  of  thinning  is  related  to 
the  result  by  Bolshev  C6]  discussed  later  in  this  section.  Lewis  and 
Shedler  C143  have  developed  an  algorithm  for  nonhomogeneous  Poisson 
point  processes  which  is  also  related. 

Ahrens  and  Dieter  C23  proposed  algorithm  PG  which  uses 
relationships  between  the  Poisson,  gamma  and  binomial  distributions  to 
generate  Poisson  variates  in  time  increasing  with  in(y>.  In  their 
computational  results,  the  execution  time  is  greater  than  for  other 
algorithms  unless  the  mean  is  quite  large.  However,  newer  algorithms 
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for  gamma  generation  (see,  e.g.,  Cheng  C83,  Schmeiser  and  Lai  C19D)  and 
binomial  generation  (see,  e.g.,  Devroye  and  Naderisamani  C10D)  make  this 
algorithm  more  competitive. 

Ahrens  and  Dieter  C2D  also  developed  a  third  algorithm  based  on 
composition.  In  the  Ahrens  and  Dieter  algorithm  PT,  a  triangular 
density  is  used  to  return  the  variate  most  of  the  time.  The  other  parts 
of  the  distribution  are  more  time  consuming  but  occur  infrequently.  The 
execution  time  increases  with 

The  acceptance/ rej ection  algorithm  is  the  basis  for  three  recent 
Poisson  generation  algorithms,  all  of  which  have  execution  times  which 
do  not  increase  (and  in  fact  decrease  slightly)  as  p  ®.  The 
acceptance/rejection  algorithm  centers  on  a  function  t(x)  which 
majorizes  f(x),  the  density  function  from  which  variates  are  to  be 

oo 

generated.  The  density  function  r(x)  =  t(x)/  f  t(y)dy  is  proportional 

•-00 

to  t(x>.  The  acceptance/rejection  algorithm  is 

1.  Generate  x  ^  r(x) . 

2.  Generate  v  *  U(0,1). 

3.  If  v  <  f(x)/t(x),  then  return  with  x  as  the  generated  variate. 
Otherwise,  go  to  Step  1,  thereby  rejecting  x. 

The  selection  of  any  function  t(x>  satisfying  t(x)  f(x)  for  all 

x  e  (-•,•)  yields  a  valid  algorithm.  Whether  the  algorithm  is  good 
depends  upon  the  speed  of  performing  Step  1,  the  difficulty  in 
evaluating  the  ratio  in  Step  3,  and  the  expected  number  of  iterations 


required  to  generate  one  variate.  Atkinson  [53  proposes  algorithm  PA 
which  uses  a  logistic  majorizing  function  and  Devroye  [93  proposes 
algorithm  JP  which  uses  a  normal  majorizing  function  for  the  body  of  the 
distribution  and  exponential  distribution  for  the  right  tail.  Algorithm 
PA  uses  tabled  values  for  x!  for  x=0,1, . . .,200.  Algorithm  ^P  uses 
preliminary  comparisons  to  avoid  calculating  x!  so  often  that  when 
evaluation  of  x!  is  required,  it  is  performed  explicitly  as 
x=x(x-1)...(3)(2).  Ahrens  and  Dieter  [33  develop  an  algorithm  based  on 
a  double  exponential  majorizing  function. 

Kronmal  and  Peterson  [133  describe  the  "acceptance/ complement" 
method,  which  is  a  composition  approach  which  requires  one  region  to  be 
generated  using  acceptance/ rej ection.  Set-up  time  can  be  reduced  by 

forcing  the  probability  of  rejection  to  be  equal  to  the  probability  of 

generating  a  variate  from  the  second  composition  region.  Ahrens  and 

Dieter  [43  develop  an  acceptance/complement  algorithm,  KPOISS,  based  on 
the  normal  distribution,  that  dominates  their  earlier  algorithm  in  [33. 

Four  approaches  which  provide  variates  which  are  approximately 
Poisson  have  been  proposed.  Atkinson  [53  includes  the  approach 
developed  in  Marsaglia  [153  and  Norman  and  Cannon  [163  which  is  based  on 
composition  and  tabling  many  values.  It  inherently  requires  PO(=x}  to 
be  truncated,  although  the  amount  of  truncation  may  be  limited  by 

increasing  the  table  size.  This  algorithm  could  be  considered  when 
memory  is  not  a  problem,  a  small  error  is  acceptable,  and  many  Poisson 
variates  are  to  be  generated  for  a  fixed  value  of  y. 

The  second  approximate  approach  is  to  use  a  normal  approximation  to 
the  distribution.  Pak  [173  discusses  the  normal  approximation  to  the 
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distribution  of  X,  (X  +  .375)^^,  and  (X  -  1/24)^^,  where  X  is  the 
Poisson  random  variable. 

The  third  approximate  approach  is  to  use  Walker's  C22D  alias 
method.  The  method  requires  truncation  of  the  right  tail  of  the 
distribution,  memory  requirements  increase  linearly  with  the  mean,  and 
set-up  time  is  substantial  for  large  values  of  the  mean.  An  alias 
algorithm  was  the  fastest  method  for  generating  Poisson  variates 
according  to  Atkinson  C5].  This  approach  could  be  made  exact  by  using  a 
composition  framework  to  obtain  the  tail  values. 

The  fourth  approximate  procedure  is  based  on  an  exact  result  by 

Bolshev  C63:  If  (X^,X^, . . .,Xn)  is  a  multinomial  random  vector  with 

parameters  y  and  p.  =  1/n  for  i=1,2,...,n  and  y  is  a  Poisson  random 

variable  with  mean  ny,  then  X,j  ,X£, . .  .,Xn  are  independent  Poisson  random 

variables  each  with  mean  u.  Tadikamalla  [21]  suggested  using  the  normal 

distribution  to  generate  y,  noting  that  the  error  can  be  made 

arbitrarily  small  by  selecting  n  large.  Despite  the  constant  execution 

time  of  generating  y  from  the  normal  distribution,  the  algorithm's 

execution  time  as  implemented  by  Tadikamalla  increases  linearly  with  u. 

However,  the  existence  of  a  multinomial  algorithm  with  execution  time 
n 

robust  to  y  =  i  X.  would  make  the  use  of  Bolshev's  result  very 
i=1 

appealing.  Note  that  Bolshev's  result  can  be  used  to  create  an  exact 
algorithm  by  generating  y  by  algorithm  PA,  IP,  KPOISS,  or  PTPE. 

3.  ALGORITHM  PTPE 

The  Poisson  random  variate  generation  algorithm  PTPE  is  developed 
in  this  section.  Generation  of  variates  is  via  acceptance/rejection. 
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based  on 

f(x)  =  u(y-M)M!/y!  -0.5  <  x  <  •  Cl) 

=  0  elsewhere, 

where  M=<p>,  y=<x+.5>,  and  <  s  >  denotes  the  integer  portipn 
of  s.  The  function  fix)  is  constructed  by  rescaling  the  Poisson 
probability  function  f^ly)  by  the  value  of  the  function  at  the  mode  M. 

This  specific  scaling  has  three  advantages: 

1.  f  IM)  =  1  for  all  p,  thereby  reducing  set-up  time. 

2.  Machine  accuracy  evaluation  of  f(y)  requires  fewer  terms  of 
Stirling's  approximation  than  does  f^ly),  because  the  errors 
in  M!  and  Y!  tend  to  cancel. 

3.  fix)  is  numerically  stable. 

Although  details  will  remain,  specification  of  the  majorizing 
function  t(x)  and  minorizing  function  blx)  defines  the  basic  structure 
of  the  algorithm  as  shown  in  Figure  A.  The  majorizing  function  is 

k^expC-x^lx^  -  x  -  .5)3  -«  <  x  £  x^-.5 

tlx)  =  (1  +  c)  -  |M  -  x  +  . 5 1 /p^  xl~.5  <  x£  xR-.5  12) 

c  expC-xRlx  +  .5  -  xR)3  x  >  xR-.5 

and  the  minorizing  function  is 

1  —  | M  —  x  +  .5|/p.  x-.5  <  x  <  x-.5  (3) 

b(x)  =  l  L  R 

0  elsewhere. 

The  constants  kL,  x^,  xR,  c,  p^,  x^,  and  xR  are  defined  as  functions  of 
p  in  the  set-up  step  of  the  algorithm.  Proposition  1  in  the  Appendix 
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addresses  the  validity  of  tlx)  as  a  majorizing  function  of  fix). 

Figure  A  about  here 

Composition  based  on  four  regions  l subdensi ties)  is  used  to 
generate  variates  from  the  density  function  proportional  to  tlx). 
Region  1 ,  which  is  the  area  under  blx),  is  triangular  with  zero 
probability  of  rejection.  Region  2  contains  the  two  parallelograms 
which  can  be  generated  as  uniform  variates.  Regions  3  and  4  are 
negative  exponential.  Py  p.,,  Py  and  p^  in  the  set-up  step  are  the 
cumulative  values  needed  to  randomly  select  the  region  to  be  used  in 
each  iteration.  The  probability  of  selecting  each  region  is 
proportional  to  its  area. 

Algorithm  PTPE  Ij,  >  1_0 ) 

Step  0.  ISet-up  constants  as  function  of  p.  Execute  whenever  the 
value  of  p  changes.) 

M  =  <p>,  p1  =  <2.195  /M  -2.2>+0.5, 

c  =  0.133  +  8.56/16.83  +  p), 

xM  =  M+0.5,  xL  =  xm  -  p y  xR  =  xm  +  Py 

a  =  Ip  -  xL)/p,  *L  =  all+a/2), 

a  =  lxR  -  p)/xR/  xR  =  all+a/2), 

P2  =  p1 11  +2c) , 

P3  =  p2  +  10. 109+8. 25/110. 86+p))/^ 

P4  =  p3  +  C/XR' 


Step  1.  iBegin  logic  to  generate  next  variate.  Generate  u  for  selecting 


the  region.  If  region  1  is  selected,  generate  triangularly 
distributed  variate  and  return.) 

Generate  u  «/>  U(0,p^),  v  v>  l)(0,1). 

If  u  >  p.,  go  to  2.  Otherwise  return  y  =  <x„-p.v+u>. 

*  n  I 

Step  2.  (Region  2.  Parallelograms.  Check  whether  Region  2  is  used. 

If  so,  generate  y  uniformly  in  CxL-.5,  xR-.5D  and  go  to  Step  5 
for  acceptance/ rej ection  comparison.) 

If  u  >  p2,  go  to  3. 

Otherwise  x  =  x^  +  (u-p^/c, 

v  =  vc  +  1  -  |M-x+0.5|/p1 . 

If  v  >  1 ,  go  to  1 . 

Otherwise  set  y  =  <x>  and  go  to  5. 

Step  3.  (Region  3,  Left  tail) 

If  u  >  p^,  go  to  4. 

Otherwise  set  y  =  <xL  +  ln(v)/XL>, 

If  y  <  0,  go  to  1 . 

Otherwise  set  v  =  v(u-p2)XL  and  go  to  5. 

Step  4.  (Region  4,  Right  tail) 

Set  y  =  <xR  -  tn(v)/XR>, 
v  =  v(u-p3)xR. 
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Step  5.  (Acceptance/Rejection  comparison) 

5.0  (Test  for  method  of  evaluating  f(y)) 

If  M  >  100  and  y  >  50,  go  to  5.2. 

5.1  (Evaluate  f(y)  via  the  recursive  relationship 

f (y)=f (y-1 )p/y.  Start  the  search  from  the  mode.) 

F  =  1.0 

If  M  <  y. 

Then  set  I  =  M  and 
Repeat 

1=1+1 
F=F  u/I 
until  I=y. 

Otherwise 

If  M  >  y. 

Then  set  I=y  and 
Repeat 

1=1+1 
F=FI/p 
unt  i  l  I=M 

End  if . 

End  i  f 

If  v  >  F,  go  to  1. 

Otherwise  return  y. 

5.2  (Squeezing,  check  the  value  of  in  v  against  upper 
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and  lower  bounds  of  in  f(y).) 

x  =  y. 

q  =  C  jrx)/x. 

Ug  =  x-u+(x+.5)q(1+q(-.5+q/3))  +  .00084. 

A  =  tn(v) . 

If  A  >  Ug,  go  to  1 . 

0  =  (x+.5)  q4/4. 

If  q  <  0,  D  =  D/(1+q) . 

If  A  <  Ug  -  D  -  .004,  return  y. 

5.3  (Perforin  final  acceptance/ rej ection  test  by  using  the 
expression  of  m  f(y)  derived  from  the  Stirling's 
formula.) 

If  A  >  C(M+.5)tn(M/u)  +  Cx+.5Hn(M/x  -  M  +  x 
+  (1./M  -  1./x)/12. 

+  C1./x3  -  1./M3)/360.D  go  to  1. 

Otherwise  return  y. 


Remark  1_. 

In  Step  1,  Region  1  is  selected.  Since  Region  1  lies  entirely 
under  f(x),  the  probability  of  rejection  is  zero.  Since  u  j>U(0,p^), 
then  u/p^  vU(0,1).  The  triangularly  distributed  variates  are  generated 
as  the  sum  of  two  independent  uniform  variates,  denoted  w  and  v.  Then 

y  -  M  +  0.5  +  (w+v-1 )p^ 

=  +  wp^  -  (1-v)pr 


Since  v  s  11(0,1),  then  (1-v)  is  also  0(0,1).  Replacing  w  by  u/p^  and 
(1-v)  by  v  yields  the  expression  used  in  Step  1. 

Remark  2_. 

In  Region  2,  x  is  uniformly  distributed  between  M-p^  and  M+p^. 
Since  u  is  uniformly  distributed  between  p^  and  Pp  in  this  region, 
w  =  (u-p^ l/Cp^-p^ )  is  1 1(0,1).  From  the  setup,  Pp-p ^  is  equal  to  2cpy 
Substituting  into  x  =  x^  +  2wp^  yields  the  expression  for  x  used  in  Step 
2.  The  expression  for  v  results  in  vs  U(b(x),  b(x)+c),  where 
b(x)  =  1  -  |M-x+.5|/p^  is  the  triangle  of  Region  1. 

Remark  3.. 

In  Step  3,  x  is  the  negative  of  a  negative  exponential  random 
variate.  The  upper  bound  of  x  is  x^  and  the  mean  is  xL  -  1/xL* 
Similarly  in  Step  4,  x  is  negative  exponentially  distributed  with  lower 
bound  Xp  and  mean  xR  +  1/xR. 

In  Region  3,  the  accept/reject  variate  v  s  U(0,t(x>),  is 

v  =  wkL  expC-xL(xL  -  x  -  .5)3  where  w  .r  11(0,1). 

The  exponential  variate  x  can  be  generated  as  x  =  xL  -  0.5  +  in(v')/xL, 
where  v'  ^U(0,1).  Then  v  =  wk^xpC-x^-anlv* ) / xL) □  =  wk^'.  Replacing 
w  by  (u-Pp) /(p^-pj)  and  (p^-pp)  by  k^/ x^  yields  the  result  in  Step  3.  A 
similar  derivation  leads  to  the  expression  used  in  Step  4. 

Remark  4^. 

In  Step  5.0,  a  test  is  made  to  select  the  method  of  evaluating 
f (y) .  The  criteria  used  here  is  based  on  both  M  and  y.  For  small 
values  of  M  and  y,  direct  calculation  using  the  recursive  formula  is 
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faster  than  evaluating  the  bounds  derived  from  the  Stirling's  formula. 
Step  5.1  is  similar  to  algorithm  P£  by  Fishman  with  f(y)  in  place  of 
fp<y),  but  requires  no  tabled  constants.  In  Step  5.2,  a  preliminary 
test  is  made  by  comparing  inlv)  against  upper  and  lower  bounds  of 
in  f<y) .  The  expressions  of  in  f (y)  and  its  bounds  are  given  in  the 
Appendix. 

Remark  5_. 

The  idea  underlying  the  setup  for  Region  3  is  to  pass  the 
majorizing  function  t^x)  through  the  point  f(xL~.5)  and  f(xL~1.5). 
This  same  approach  is  used  in  Region  4  using  flxR-0.5)  and  flxR+0.5). 
The  result  is  t^<x)  in  Figure  B. 

Figure  B  about  here 

This  exact  set-up  requires  six  logarithms  and  two  exponential 
operations.  These  operations  are  slow  and  can  be  avoided.  The  setup  in 
PTPE  uses  the  majorizing  function  tlx)  as  shown  in  Figure  B,  which  does 
not  require  higher  order  operations.  The  use  of  tlx)  increases  the 
probability  of  rejection  slightly,  but  the  gain  in  efficiency  by 
avoiding  higher  order  operations  in  the  setup  is  significant  in  the 
cases  where  the  value  of  mean  p  changes  often.  That  these  exponential 
tails  majorize  fix)  is  proved  in  the  Appendix. 

Remark  6. 

The  expected  number  of  U10,1)  values  required  to  generate  a  Poisson 

_  M 

variate  is  2lp^)le  11  w  /M!),  where  M  is  the  integer  portion  of  w  and 
P 4  *  J tlx)dx,  as  defined  in  Step  0. 


The  derivation  is 


straightforward.  The  expected  number  of  iterations  is 

t ( x) dx/  /  f(x)dx=  p J  (M! /e~ p  UM)  /fp(x) 

”  oo  ™  oo 

=  P4(e"‘1  yM/M!). 

Multiplying  by  the  two  U(0,1)  values  per  iteration  yields  the  result. 
Remark  7_. 

All  four  fundamental  concepts  are  included  in  PTPE.  The  overall 
structure  of  PTPE  is  acceptance/ rej ection.  The  inverse  transformation 
is  used  to  select  the  region,  to  generate  uniformly  distributed  variates 
for  Region  2,  and  to  generate  exponentially  distributed  variates  for 
Regions  3  and  4.  The  use  of  four  regions  is  composition.  The  special 
property  that  the  sum  of  two  independent  U(a,b)  random  variables  has  a 
triangular  distribution  is  used  in  Region  1. 

4.  COMPUTATIONAL  EXPERIENCE 

The  four  algorithms  for  which  execution  time  approaches  a  constant 
as  u  ■»  ®,  PA,  IP,  KPOISS,  and  PTPE,  are  compared  here  in  terms  of  setup 
times,  marginal  execution  times,  lines  of  code,  and  memory  requirements. 
The  Ahrens  and  Dieter  algorithm  in  [33  is  dominated  by  KPOISS  and  not 
discussed  here.  All  four  algorithms  were  implemented  in  FORTRAN  using 
the  MNF  compiler  on  Purdue  University's  CDC  6500  computer.  The  uniform 
(0,1)  variates  were  generated  using  RANF,  which  is  intrinsic  in  the  MNF 


as 

/ 
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For  each  combination  of  y  and  algorithm,  four  replications  of  3000 
variates  were  timed.  The  execution  times  shown  in  Table  1  are  the 
averages  of  the  replication  averages  and  are  accurate  to  within  one  unit 
in  the  last  decimal  place.  The  accuracy  may  also  be  assessed  by 
comparing  the  last  four  lines  in  the  table,  for  which  most  of  the 
differences  in  times  are  due  to  random  variation  rather  than  to  changes 
in  distribution  shape. 

The  marginal  execution  times,  shown  under  the  heading  "Fixed  Mean" 

in  Table  1,  favor  PTPE.  The  execution  times  for  setting  up  the 

algorithm  and  generating  one  variate,  shown  under  the  heading 

-9 

"Incremented  Means"  in  Table  1,  were  obtained  by  incrementing  y  by  10 
with  each  variate  generated.  Because  KPOISS  requires  little  more  than  a 
square  root  calculation  to  set  up,  it  is  competitive  with  PTPE  when  the 
mean  changes  with  each  variate  generated. 

Since  IP  and  KPOISS  require  normal  variates,  their  times  are 
sensitive  to  the  normal  variate  generation  algorithm  used.  We  used 
algorithm  KJR  (see  Kinderman  and  Ramage  Cl  2D  >  which  is  the  fastest 
FORTRAN  level  algorithm  available.  For  those  who  have  a  faster 

assembler  language  normal  generator  available,  the  times  for  KPOISS  and 
IP  would  be  less.  Of  course,  all  four  algorithms  would  be  faster  if 
coded  in  assembler  language.  Another  comment  concerns  PA.  The 
approximation  given  by  Atkinson  C  5  3  for  the  constant  c  is 

c  =  .767  -  3.36/ g,  which  is  inaccurate  when  y  <  30.  The  poor 
approximation  causes  the  relatively  large  execution  times  of  PA  for 

u* 


small  values  of 


17 


Table  1.  Comparison  of  Algorithms 


Fixed  Mean 

Incremented 

Mean 

PTPE 

KPOISS3 

IP3 

PA 

PTPE 

KP0ISS 

IP 

PA 

iob 

.33 

.38 

.66 

1.41 

.35 

.45 

1.34 

1.78 

25 

.29 

.37 

.62 

1.03 

.50 

.45 

1.31 

1.41 

100 

.24 

.36 

.58 

.90 

.48 

.44 

1.27 

1.27 

250 

.22 

.35 

.57 

.89 

.44 

.43 

1.26 

1.26 

1000 

.20 

.34 

.55 

.86 

.42 

.42 

1.25 

1.24 

10,000 

.20 

.34 

.54 

.87 

.42 

.42 

1.23 

1.24 

1,000,000 

.20 

.34 

.54 

.87 

.41 

.41 

1.23 

1.23 

Memory 
Requi re¬ 
men  ts 

279 

309 

282 

146 

Lines  of 
Code 

59 

64 

64 

17 

a 

Using 

KR  for  the 

normal 

random 

variates. 

b 

Times 

for  KP0ISS 

are  for 

w  =  10 

+  G . 
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The  execution  times  were  compared  using  a  slower  U(0,1)  generator. 
All  times  in  Table  1  increased  by  about  .1  except  for  which  had 
increases  of  about  .2. 

The  execution  times  were  also  compared  using  the  FTN  compiler.  For 
large  values  of  p,  KPOISS  required  only  56%  more  time  than  PTPE 
(compared  to  75%  under  MNF)  for  fixed  means.  For  variable  means  KPOISS 
was  9%  faster  than  PTPE  (compared  to  0%  under  MNF). 

Note  that  several  exact  algorithms  are  faster  than  the  four 
algorithms  compared  here  for  small  values  of  the  mean  (approximately 
M  <  50). 

While  the  number  of  lines  of  FORTRAN  code  is  only  a  crude  measure 
of  the  goodness  of  an  algorithm,  it  can  be  important  both  in  terms  of 
the  effort  to  implement  the  algorithm  and  to  verify  that  the  algorithm 
is  working  properly.  PA,  PTPE,  KPOISS,  and  ^P  required  17,  59,  64  and 
64  lines  of  code,  respectively.  This  does  not  include  the  nine  lines 
for  the  routine  used  to  evaluate  *n(x!)  needed  by  PA  nor  the  58  lines  of 
the  KR  normal  variate  generator  used  here  by  ^P  and  KPOISS .  Algorithms 
PA,  PTPE,  IP,  and  KPOISS  require  146,  279,  282,  and  309  words  of  memory, 
respectively,  again  not  including  required  support  routines. 
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APPENDIX:  PROPERTIES  OF  b(x)  AND  tlx) 

Four  inequalities  used  in  algorithm  PTPE  are  discussed  here. 
Proposition  1  considers  b(x)  _<  f(x)  <  t(x>,  which  is  necessary  for  the 
acceptance/ rej ection  parts  of  PTPE.  In  addition,  in  PTPE  f(x)  is 
squeezed  by  upper  and  lower  bounds  which  are  pi oved  valid  in 
Propositions  2  and  3,  respectively. 

Results  1-3,  stated  below  without  proof,  are  necessary  for  the 
proofs  of  Propositions  1,  2,  and  3.  All  follow  from  the  Taylor  series 
expansion  of  the  logarithm  (see  e.g.,  Abramowitz  and  Stegun  Cl]). 

Result  1_.  If  a  <  b,  then  ln(b/a)  >.  d  +  q2/2,  where  q  =  (b-a)/b. 

Result  2.  For  all  a  >  0  and  b  >  0,  *n(b/a>  <  q  -  q^/2  +  q^/3,  where 

q  =  (b-a)/a. 

Result  3^.  For  all  a  >  0  and  b  >  0,  *n(b/a>  *1  q  -  q^/2  +  q^/3  -  A q^/4, 

where  q  =  (b-a)/a  and  A  =  1  if  a  £  b  and  A  =  (1+q)  **  if  a>b. 

Lemma  1  is  used  in  the  proof  of  Proposition  1. 

Lemma  1.  For  all  u  >  0, 

£  f^lx)  £  fM<x>  if  x=0,1,2, . . .,M 

and 

f*+1-e<x>  >  fp<x>  >  f^Cx)  if  x=H,M+1 , . . ., 
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where  M  =  <w>  and  fy(x)  =  ux  MM!/x!. 

Proof.  The  ratio  f  (x)/f„(x)  =  (u/M)x  Since  M  <  u,  the  right  side 

-  U  n  — 

inequalities  follow.  Similarly,  the  left  side  inequalities  follow  from 
f^.  (x)/f  (x)  =  C(M+1-e)/u)x"M.  || 

H+l—t  y 

Proposition  1_.  For  w  10  and  xe  (-«,<»),  b(x)  <  f  (x)  <  t(x) ,  where  t(x) 
and  b(x)  are  defined  in  Equations  (2)  and  (3),  f(x)  is  defined  in 
Equation  (1),  and  specific  constants  are  defined  in  Step  0  of  algorithm 
PTPE . 


Proof .  The  proof  is  trivial  for  xc  (-*,-.5),  since  f(x)  =  0. 
Consider  xe  (-.5,xL  -.5).  Since  xL  >  0, 

( x^/x^)  (x^/lx^-1 )) .  ,.(xL/<  x+1.5  >)  ^  1 


which  implies 


Then 


x  -<  x+.5  > 

xL  21  xl!/<  x+.5  >!  • 


x. -<  x+.5  >  ,  . 

xlL  >  (xL/v)<x  *5>“<  x  '5  xL ! /<  x+.5  > ! 


since 

Direct  algebra  yields 


xL  <u  and  (x+.5)  <  x+.5  >. 


which  implies 


x  -<x+.5)  <  x+.5  >-M-(x.  -M)  +  (x  -(x+.5)> 

XL  >  »  xl!/<  x+.5  > ! 


x.  -M  x.  -(x+.5)  x  c  N  M 

Cu  L  M!/xl!)(xl/m)  L  >  m  x  M  M! /<  x+.5  >! 


which  implies 


f(xL)  expC(xL-(x+.5>) tn(xL/M)D  >  f(x). 


(A-1) 


Applying  Result  1  to  in(x^/p)  =  -^ntP/x^)  yields 

f(xL)  exp[- Al(x^-(x+.5))3  >  f(x),  ( A— 2 ) 

where  A  =  a.  +  a2/2  with  a.  =  C p—  x ,  ) / y . 

L  L  L  L  L 

The  majorizing  function  used  in  the  algorithm,  valid  for  p  £  10,  is 

k^  exp[- A^(x^-(x+. 5) >D  £  fix),  (A-3) 

where  k^  =  .109  +  8.25/(10.86  +  p> .  Inequality  (A-3)  requires 

f(xL>  £  kL  for  all  p  £  10.  Since  k^  £  M,  Lemma  1  implies  that  only 

integer  values  of  u  need  be  considered.  The  inequality  was  numerically 

verified  for  p  =  10,11 ,.. .,10000.  The  proof  that  f(x^)  £  k^  for 

v  e  C10000, ”)  is  based  on  showing  ffx^)  £  z^(xL>  £  z^lx^)  £  k^,  where 

Zi(xL)  =  exp[-(xL~p)2/2p]  and  z^(xL)  =  exp[~(2.195  J~ v  -3.2)2/(2p)]. 

The  left  inequality  is  from  the  normal  majorizing  function  used  by 

Ahrens  and  Dieter  [4D  for  all  x  £  <  p-1.1484  >.  The  center  inequality 

follows  from  -(xl-p)  =  -<  2.195  /p  -2.2  >2  £  -(2.195  ft  -3.2)2.  The 

right  inequality  follows  from  z^dOOOO)  =  .0964  <  min  kL  =  .109  and  that 

p 

for  all  p  e  [10000,°°),  z.,(xL)  3  decreasing  function  of  p.  That 

-3/2  -2 

z.,(xL)  decreases  fotlows  from  d  tn  z^fx^/dp  =  -3.512P  +  5.12p 

which  is  negative  for  all  p  >  2.1254. 

Similar  logic  for  x  e  (xR  -.5,®)  leads  to 

c  exp[-AR(x+.5-xR)D  £  f(x),  (A-4) 

where  c  =  .133+8. 56/(6. 83+p)  for  alt  p  £  10. 

Now  consider  x  e  [xL-.5,xR-.S3,  for  which  b(x)  £  f(x)  £  t(x)  must 
be  satisfied,  where  b(x)=  1  -  |M-x+.5|/p^  and 
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t(x)  =  (1+c)  -  |M-x+.5 | /p^ .  Again,  y  e  CIO, 10000]  was  checked 

numerically;  using  y  =  M+1-e  when  x  £  M  and  y  =  M  when  x  £  M  for 

b(x)  £  f(x)  and  y  =  M  when  x  £  M  and  y  =  M+1-e  when  x  £  M  for 

fCx)  £  tlx),  as  indicated  by  Lemma  1.  For  y  £  10,000,  consideration  of 
limiting  values  and  the  asymptotic  value  of  .133  for  c  indicates  the 
inequality  is  satisfied.  || 

Lemmas  2,  3,  and  4  are  needed  for  the  proofs  of  Propositions  2  and 
3,  which  are  upper  and  lower  bounds  on  f(x),  respectively. 

Lemma  2.  For  M  =  <  y  >,  (M+. 5) *n(M/ y)  £  M  -  y. 

Proof.  Substituting  x=M/y  into  the  well-known  inequality  fcn  x  £  x-1 
and  multiplying  by  M+.5  yields  (M+. 5) in(M/ y)  £  CM/ y) (M-y)+ (M-y)/(2 y) . 
Since  (M-y)/(2  y)  £  0  and  0  £  M/  y  £  1 ,  the  result  is  obtained .  | | 

Lemma  3.  For  y  £  y*  and  M  =  <  y  >, 

M  -  y  +  g(y*)  £  (M+. 5) *n(M/ y) , 

where  g(y*)  =  «  y*  >  +  . 5) C ln( y*/( y*  +  .5))]  +  .5. 

Proof.  The  proof  shows  that  g(y*)  minimizes 

g(y)  =  Min  C (M+. 5) ln(M/ y>- (M- y)] . 

* 

v  <  y 
M=  <y> 

First  consider  Min  g(y).  betting  dg(y)/dy  =  0  and  checking  that 
2  2  M  £  y  <  M+1 

d  g(y)/dy  >  0  yields  y  =  M+.5.  The  problem  is  now  to  find  the  value  of 
M  which  minimizes  (M+.5) tn(M/ (M+. 5) )+. 5  subject  to  M  £  <y*>.  Since  the 
function  increases  with  M,  as  can  be  seen  graphically  or  by  evaluating 
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derivatives,  the  optimal  value  is  M  =  <  u*  >.  || 

Lemma  4.  Consider 

6<M.y)  =  (M_1-y_1)/12  -  (M~3-y“3 ) /360  +  (M_5-y ~5) /I  260, 
«L(M*,y*)=  -(12y*)-1  -  (360M*3)-1  -  C1260y*5)-1, 

and 

«u(M*,y*)=  (1 2M*)-1  +  (360y*?)_1  +  (1260M*5)”1. 

IF  M  >  M*  and  y  >  y*,  then  1  6u<M*,y*). 

Proof .  The  lower  and  upper  bounds  are  obtained  directly  by  minimizing 
term  by  term  for  5L(M*,y*>  and  maximizing  term 
by  term  for  ^(M^y*).  || 

Proposition  2.  Consider 

=  y-u+(y+.5)q(1+q(-.5+q/3))+ 6  (M*,y*) 

where  q  =  (  m- y)/y.  If  M  j>  M*  and  y  y* ,  then  l)^  >.  in  f(y). 

Proof.  The  proof  algebraically  simplifies  Jn  f(y),  which  is  evaluated 
using  Stirling's  Formula.  Further  simplification  results  from 
inequalities  on  relatively  insignificant  terms. 
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in  f  (y)  =  (y-M)  In  y  t  in  M!  -  in  y! 

=  (y-M)  in  u 

+  C(M+.5)ln  M  -  M  +  In  +  (12M)"1 

-  (360M3)'1  +  (1260M5)"1  +  o(M)“7] 

-  C(y+.5)in  y  -  y  ♦  In  f2*  +  (1 2y)_1 

-  (360y3)-1  +  (1260y5)-1  +  o(y)“7] 

=  C(y+.5)  -  (M+.5)Din  u  +  (M+.5)in  M  -  (y+.5)£n  y 

-  M  +  y  +  6(M,y) 

=  (M+.5)  in(M/p)  +  (y+.5)in(p/y)  -  M  +  y  +  6(M,y)  (A-5) 

where 

«(M,y)  =  (M_1-y-1)/12  -  (M~3-y"3)/360  +  (M-5-y~5) /1 260  +  o(CM-y)-7). 

Applying  Lemma  2  to  the  first  term.  Result  2  to  in( y/y)  and  Lemma  4  to 
A(M,y)  in  Equation  (A-5)  yields  the  result. 

Proposition  3.  For  M  £  M*  and  y  >  y* , 

-  D  +  g(p*)  -  ^j(M*,y*)  -  d^(M*,y*)  £  tn  f(y), 

where  D  =  (y+.5)q^A/4,  A  =  1  if  q  >  0  and  A  =  (1+q)-^  if  q  <  0,  and 
q  =  ( ir-y) /y. 


27 


Proof.  Ub  -  0  +  g(u*>  -  d^CM^y*)  -  ^(M^y*) 

=  Cy  -  m  +  (y+.5)q(1+q(-.5+q/3))  +  ^(M*,y*)]  +  g(u*> 

-  C(y+.5)q4  A/A3  -  «,j(M*,y*>  +  ^(M^y*) 

=  y  -  u  +  (y+.5)Cq-q^/2+q^/3-q^ A/A:  +  g(y*)  +  fi^CM^y*).  (A-6) 

From  Lemma  3,  g(y*)  <  (M+. 5)  4n(M/  u)  -  M  +  vi;  from  Lemma  A, 
6^**'^*^  5.  «™,y) ;  and  from  Result  3  applied  to  tn(u/y>.  Equation  (A-6) 
is  less  than  in  f (y) . | | 
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Computer  codes  used  to  obtain  the  conputational  results 
of  section  4,  "Poisson  Random  Variate  Generation"  by 
Bruce  Schmeiser  and  Voratas  Kachitvichyanukul. 


THIS  IS  THE  MAIN  PROGRAM  TO  TEST  CARIOUS  METHODS 
G~  GENERATING  P0ISS0M  RANDUI1  UARIATES 

ERUCE  SCHMEISER  AMD  UORATAS  KACHITUICHYANUKUL 
SCHOOL  OF  INDUSTRIAL  ENGINEERING 
PURDUE  UNIUERSITY,  AFRIL.  1330 


DIMENSION  NAM£(5)»XXMU( 13) 

DATA  NAME/’ CUNY’ . ’ PTPE’ , ’ KPOISS* . ’ IP’ . ’ PA’ / 

DATA  XXMU/10. ,25. , 50. , 50.5, 100. .230.. 500., 1000., 3000., 

1  5000. ,10000. ,50000. .1000000./ 

N=3000 
ISEED=0 
kRITECG, 1000) 

1000  FORMAT ( 1H1 ) 

DO  400  L=l» 13 
XMU=XXMU(L) 

11=0 

Ti1EAN=XXMU(L) 

TUAR=XXMU(L) 

STE=S0RT  C  TUAR/N ) 

URIT£(6»3000)  N. TMEAN. TUAR. STE 

3000  F0RMATC3X,  ’  »«s*s»ssiicesH»*u»»»#»»*3***»*»**»**»*****************’ 

1  /3X, SAMPLE  SIZE  =  ’15.’  :::::'/ 

2  /11X, ’TIME  MEAN  UARIANCE  STD  ERROR’/ 

3  ’  TRUE  ’.3F15.3) 

DO  100  1=1,5 

SUMT=0. 0 

DO  150  J=l, 4 

SUM=0.Q 

SUM2=0.0 

CALL  SECOND (Tl) 

C 

DO  300  K=1.N 

CO  TO  (1,2. 3. 4,5). I 

1  CONTINUE 
CO  TO  200 

2  CALL  PTPECXMU. ISEED, II) 

GO  TO  200 

3  CALL  KPOISSCXMU. ISEED. II) 

GO  TO  200 

4  CALL  IP(XMU. ISEED. II) 

GO  TO  200 

5  CALL  PA(XMU. ISEED. II) 

200  SUM=SUM+I I 

SUM2=SUM2+ I I » I I 
300  CONTINUE 

C 

CALL  SECOND (T2) 

TIME=1000.»(T2-T1)/N 

SUMT=SUMT+TIME 

AUCT-SUMT/J 

XMEAN*SUM/N 

UAR*SUM2/N-XMEAN*XMEAN 
URITE(G,2000 )  NAMEC I ) . TIME, AUGT.XMEAN, UAR 
2000  FORMAT (IX , A5 , 2F 6 . 3, F 15 . 3, F 15 . 3 ) 

150  CONTINUE 

100  CONTINUE 
400  CONTINUE 
C 

STOP 

END 
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SUBROUTINE  PTFEQ\MU. ISEED.  JX) 

C 

C  POISSON  RANDOM  UARIATE  GENERATOR 

C  M!U  :  MEAN  CXirj  .  G£.  10) 

C  ISEED  :  RANDOM  NUMBER  SEED 

C  JX  :  RANDOMLY  GENERATED  OBSERVATION 

C 

C  ERUCE  U.  SCHMEISER  AND  UORATAS  K  ACH I TUI CHYANUKUL 

C  PURDUE  UNIUERSITY.  SEPTEMBER  1930. 

C  REUISED  JULY.  1981 

C  METHOD  :  ACCEPTANCE-REJECTION  UIA  FOUR  REGION  COMPOSITION 
C  AUXILIARY  REQUIRED  SUBPROGRAM  : 

C  UNIFORM  (0.1)  RANDOM  NUMBER  GENERATOR 

C 

DATA  YMU/-1 ./ 

IF(XMU.EO.YMU)  GO  TO  B 
C 

C»****SETUP  (EXECUTE  ONLY  WHEN  XMU  CHANGES) 

C 

YMU=XMU 

M=YMU 

FM=M 

Pl=INT(2.19S*S0RT(FM)-2.2)+0.5 
C= .  1 33+8 . 5S/  ( S .  83+YMU ) 

XM=M+0.5 

XL=XM-P1 

XR=XM+P1 

AL=(YMU-XL)/'YMU 

XLL=AL»(l.+.5-*AL) 

AL=(XR-YMU)/XR 

XLR=AL*(1.+.5*AL) 

P2=P1*(1.+C+C) 

P3=P2+ ( 0 . 1 09 +8 . 25/  ( 1 0 . 86+YMU )  )  /XLL 
P4=P3+C/XLR 
C 

C«***«GENERATE  UARIATE 
C 

2  U=RANF ( ISEED ) *P4 
U=RANF( ISEED) 

C 

C  TRIANGULAR  REGION 

C 

IF(U.GT.Pl)  GO  TO  3 
IX=XM-P1*U+U 
GO  TO  14 
C 

C  PARALLELOGRAM  REGION 

C 

3  IF(U.GT.P2)  GO  TO  4 
X=XL+(U-P1)/C 
U=U>C+l.-ABS(FM-X+0.5)/Pl 
IF(U.GT . 1 . )  GO  TO  2 
IX=X 
GO  TO  G 

LEFT  TAIL 


C 

C 

C 


4  IF(U.GT.P3)  GO  TO  5 
IX=XL+ALOG(U)/XLL 
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IFCIX.LT. 0)  GO  TO  2 
U=U-(U-P2)*XLL 
GO  TO  G 

RIGHT  TAIL 

5  IX=XR-ALOGCU)/XLR 
U=U-CU-P3)*XLR 


ae*ea ACCEPT ANCE-RE JECT ION  TEST.  COMPARE  U  TO 
THE  SCALED  POISSON  MASS  FUNCTION 

G  IFCM.GE.100.fiND.IX.GT.50)  GO  TO  12 

EXPLICIT  EUALUATION 

F=1 . 0 

IFCI1-IX)  7,11,9 

7  MP=M+1 

CO  8  I=MP» IX 

8  F=F-Yi1U/I 
GO  TO  11 

9  IX1=IX+1 

DO  10  1=1X1, M 

10  F=F*I/YMU 

11  IFCU-F)  14,14,2 

SQUEEZING  USING  UPPER  AND  LOWER  BOUNDS  ON  ALOGCFCX)) 

12  X=IX 
0=CYMU-X)/X 

UB=X-YMU+  C  X+ . 5 ) «Q* ( l . +Q* ( - . 5+ . 33333333333*0 ))+. 00084 
ALU=ALOGCU) 

IFC ALU.GT.UB)  GO  TO  2 
D=CX+O.S)*.25*(0«Q)»*2 
IFCO.LT.O.)  D=D/'C1.+Q) 

IFCALU.LT .UB-D-. 004)  GO  TO  14 

STIRLING’S  FORMULA  TO  MACHINE  ACCURACY  FDR 
THE  FINAL  ACCEPTANCE^REJECTION  TEST 

I F  C  ALU . GT . C  FM+ . 5 ) «AL0G ( FM/ YMU )  +  C  X+ . 5 ) *ALOG  C  YMU/'X ) -FM+X 

1  +(1./'FM-1./XV12.  +  .0027777777778/CX*X*X) 

2  -.00277777778/(FM*FM»FM))  GO  TO  2 
14  JX=IX 

RETURN 

END 
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SUBROUTINE  KPOISSCA. IR.KPOIS) 

J.  H.  AHRENS  AND  U.  DIETER 

COMPUTER  GENERATION  OF  POISSON  DEUIATES  FROM 
MODIFIED  NORMAL  DISTRIBUTIONS 

DIMENSION  FACT(10)*PP(35) 

DATA  AA.AAA.A0.A1.A2.A3.A4,A5,A6.A7  /0.,0.,-. 5, .33333333, 

1  -.2500068, .2000118,-. 1661269. .1421878,-. 1384794,. 125006/ 

DATA  FACT  /l. . 1. ,2. ,6. ,24. , 120. ,  720. ,5040. ,40320. ,362880./ 
IF(A.EQ.AA)  GO  TO  1 
IFCA.LT. 10.0)  GO  TO  12 
AA=A 

S=SQRT(A) 

D=6.0*A»A 
L=INT(A-1. 1484) 

1  CALL  NORMAL! IR.TT) 

IF(G.LT.O.O)  GO  TO  2 
KPOIS=INT!G) 

IF(KPOIS.GE.L)  RETURN 
FK-FLOAT ( KPO IS) 

AK=A-FK 

U-RANF1IR) 

IF(D*U.CE.AK*AK«AK)  RETURN 

2  IF(A.EQ.AAA)  GO  TO  3 
AAA-A 

OMEGA-. 3389423/S 

B1-.4166667E-1/A 

B2=.3»B1»B1 

C3-.1 42857 1*B1*B2 

C2=B2-15.*C3 

C1=B1-6.*B2+45.*C3 

C0=1.-B1+3.»B2-15.*C3 

C=. 1069/A 

3  IF(G.LT.O.O)  GO  TO  5 
KFLAG=0 

GO  TO  7 

4  IF ( FY-U*FY. LE . PY»EXP ( PX-FX ) )  RETURN 

5  E--ALOG ! RANF  !  IR  )  ) 

U=RANF(IR) 

U-U-HJ-1.0 

T-1.8+SIGNIE.U) 

IF(T.LE. -0.6744)  GO  TO  5 

KPOIS-INT!A+S*T) 

FK -FLOAT ( KPO IS ) 

AK-A-FK 
KFLAG-1 
CO  TO  7 

6  IF!C*ABSCU).GT.PY*EXPCPX+E)-FY*EXP(FX+E))  CO  TO  5 
RETURN 

7  IFOCPOIS.GE. 10)  CO  TO  8 
PX— A 

PY*A**KPCI S/FACT (KPOIS+1 ) 

CO  TO  11 

8  DEL-.8333333E-1/FK 
DEL-DEL-4 . 8*KL*DEL*DEL 

U^AK/FK 
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IF(ABS(U) .LE.0.25)  GO  TO  9 
PX=FK * ALOG ( 1 . 0+U  > -AK-DEL 

9  px=FK  »U*U» (((((( ( A7*U+A6 ) *U+ AS ) *U+ A4 ) *U+A3 ) *U+A2 ) *U+A 1 ) *U+A0 ) -DEL 

10  PV*. 3989423/S0RT  C FK ) 

11  X=(0.5-AKVS 
XX=X«X 
FX=-0.5*XX 

FY=OMEGA* ( ( (C3*XX+C2)*XX+C1)*XX+C0) 

IF(KFLAG)  4»4>G 

12  AA=0.0 

IF(A.EQ.AAA)  GO  TO  13 

M=MAX0( 1« INT(A) ) 

L=0 

P=EXPC-A) 

Q=P 

P0=P 

13  U=RANF(IR) 

KPOIS=0 

IF(U.LE.PO)  RETURN 
IF(L.EQ.O)  GO  TO  IS 
J=1 

IFCU.GT. 0.458)  J=MIN0(L,M> 

DO  14  KPOIS=J»L 

IF(  U.LE.PP(KPOIS) )  RETURN 

14  CONTINUE 
IF(L.EQ.3S)  GO  TO  13 

15  L=L+1 

DO  IB  KP0IS=L.35 
P=P* A/FLOAT ( KPO I S ) 

Q=0+P 

PP(KPOIS)=Q 
IF(U.LE.Q)  GO  TO  17 

16  CONTINUE 
L=35 

GO  TO  13 

17  L=KP01S 
RETURN 
END 
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SUBROUTINE  IP(SL. ISEED. IX) 

THE  SUBPROGRAM  IP  PRODUCES  RANDOM  POISSON  UARIATES  UITH 
PARAMETER  SL>0.  A  REJECTION  METHOD  WITH  SQUEEZING  IS  USED. 

BY  LUC  DEUROYE.  MCGILL  UNIUERSITY.  CANADA  1980 

AUXILIARY  SUBPROGRAMS  REQUIRED  : 

UNIFORM  (0.1)  RANDOM  NUMBER  GENERATOR 
STANDARD  NORMAL  RANDOM  NUMBER  GENERATOR 
EXPONENTIAL  RANDOM  NUMBER  GENERATOR 

DATA  L/Ox, SLM. D. D2.  D3.STDEU. PTAIL. PBODY. CON. RL, RI . TWO. EL/12*0./ 

SETUP 

IF(SL.EQ.SLN)  GO  TO  10 
L=INT(SL) 

RL=FLOAT(L) 

RI=1./RL 
EL=EXP(RL-SL) 

TWO=RL+RL 

D=S0RT ( RL*ALOG (1 . + 1 0 . 1 8593*RL ) ) 

D2=D+TW0 
D3=D2/D 

STDEU=SQRT(0.5*D2) 

PTAIL=D3«EXPC-(D+1.)/D3) 

CON=0.25^TWO 

PBODY=EXP(CON)*SQRT(3. 141 593*D2 ) 

SUM=PTAIL+PBQDY+1 . 0 
PBODY=l./SUM 
PTAIL=PBODY+PT  AIL/SUM 
SLM=SL 
C 

10  IFCL.EQ.O)  GO  TO  99 
1  A=0. 

IX=0 

U=RANF ( 0 ) 

IFCU.LT. PTAIL)  GO  TO  50 
CALL  NORMALC ISEED. R) 

X=R*STDEU-0.5 

IF(X.GT.D.OR.X.LT.-RL)  GO  TO  1 

IF(X.GT.O)  GO  TO  18 

A=»1.0 

X=X-2. 

18  IX=INT(X+1.) 

Y=FL0AT(IX) 

U*ALOG ( RANF ( I SEED ) ) -0 . 5*R**2+C0N 

19  T=Y*(Y*1.VTW0 

IFCU.LT.-T.AND.A.EO.O.)  GO  TO  100 
QR=T»(-1.+(Y+Y+1. )*0.1BGSBG7*RI) 

OA=QR-T**2*0 . 3333333/ ( RL+ ( Y+ 1 .) *A ) 

IF(U.LT.QA)  GO  TO  100 
IF(U.GT.QR)  GO  TO  l 
RM=RI*(1.-2*A) 

K=-IX-1 

IF( IX.GT.O)  K-IX 
PD=1 .0 
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S=0. 

DO  20  J=1.K 
S=S+RM 

20  PD=PD»C 1 .  +S) 

IFCU.LT. C2»A-1. )*ALQGCPD) )  GO  TO  100 
GO  TO  1 

50  IF(U.LT.PBODY)  GO  TO  100 
X=D-ALOGCRANFC ISEED) )*D3 
IX=INT CX+1 . ) 

Y=FLOATC IX) 

U=AL0G C RANF ( I SEED ) ) - ( X+ 1 . ) / D3 
IFCU.GT.-Y»CY+1. )/CTU0+Y))  GO  TO  1 
GO  TO  19 
99  IX=0 
100  IX=IX+L 

PD=RANF ( I SEED) 

110  IFCPD.LT. EL)  RETURN 
IX=IX+1 

PD=PD*RANF ( ISEED ) 

GO  TO  110 
END 


SUBROUTINE  PACXMU. ISEED. IX) 

GENERATE  THE  POISSON  RANDOM  UARIATE  IX  UITH  MEAN  XMU 
USING  ATKINSON'S  ALGORITHM  (XMU.GE.10) 

APPLIED  STATISTICS.  28.  1C  1979).  29-39 
ALOGCIX  FACTORIAL)  UIA  STIRLING'S  APPROXIMATION 

DATA  SAUE/-1./ 

IFCXMU.EO.SAUE)  GO  TO  100 
SAUE=XMU 

B=  l .  8 1 37993S42/SORT  C  XMU ) 

A=B*XMU 

C= ALOG  C  0 . 7G7-3 . 36/XMU ) -XMU-ALOG  C  B ) 

AM=ALOGCXMU) 

100  U=RANFC ISEED) 

X=CA-AL0GCC1.-U)/U))/B 
IFCX.LE.-0.5)  GO  TO  100 
IX=X+0.5 
U=RANFC ISEED) 

Ds 1 • +EXP ( ft— B*X ) 

IFC A-B»X+ALOGCU/CD*D) ) .LE.C+IX*AM-XLFACC IX) )  RETURN 

GO  TO  100 

END 


FUNCTION  XLFACCI) 

FUNCTION  TO  EVALUATE  LOG  OF  THE  FACTORIAL  I 
I  .GT.  7  USE  STERLING’S  APPROXIMATION 
I  .LE.  7  USE  TABLE  LOOKUP 

DIMENSION  ALC8) 

DATA  AL/0. . 0. .0.6931471806. 1.791759469. 3. 1 78053830. 
I  4.787491743.6.579251212.8.525161361/' 

IFCI.GT.7)  GO  TO  100 
XLFAC=ALCI+1) 

RETURN 

100  XLFAC=C 1+0.5) *ALOG C FLOAT ( I ) )-I +0 . 08333333333333/1 
1  -0.00277777777777/C I*I»I )  ♦  0.9189385332 

RETURN 
END 
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SUBROUTINE  NORMAL! ISEED. X) 

GENERATION  OF  ONE  NORMALCO.l)  UARIATE  USING 
THE  ALGORITHM  GIUEN  BY  KINDERMAN  AND  RAMAGE 
IN  THE  JOURNAL  OF  THE  AMERICAN  STATISTICAL  ASSOCIATION  12/76 

CODED  BY  PETER  BONNER  AND  MODIFIED  BY  BRUCE  SCHMEISER 
MARCH  1977  AND  JUNE  1977  RESPECTIUELY 


DATA  TAIL/2.21S0358G71GG471/ 

UU=RANF ( I SEED) 

IFCUU.GE. .884070402298758)  GO  TO  2 

RETURN  TRIANGULAR  UARIATE  88  PERCENT  OF  THE  TIME 
Y=RANFC ISEED) 

X=TAIL*C 1.131 131635444180*UU+Y-1 .0) 

RETURN 

2  IFCUU.LT.. 973310954173898)  GO  TO  4 
TAIL  COMPUTATION 

3  U=RANF ( ISEED) 

U=RANF  C I SEED) 

T1=TAIL*TAIL'2.0 
T=Tl-ALOG(U) 

IFCU*U*T.GT.T1)  GO  TO  3 
X=SQRTC2.0»T) 

IFCUU.GE. .98GG5547708S949)  X=-X 
RETURN 

4  IFCUU.LT.. 9587208247904G3)  C.0  TO  G 
FIRST  NEARLY  LINEAR  DENSITY 

5  U=RANFCISEED) 

U=RANFC ISEED) 

Z=U-U 

LET  U=  MAXCU.U)  AND  LET  H=MINCU,H) 

IFCU.GT.U)  GO  TO  100 
TEMP=U 
U=U 
U=TEMP 

100  T=TAIL-.G308348019219G0»U 

IFCU.LE. .75559 153 IGG7G0 1 )  GO  TO  9 

DIFF=EXPC-T*T*.5)/2.50GG28274G3100-. 1800251S10G8563* 
*  C2.21G0358G71GG471-ABSCT)) 

IFCABSCZ)*. 034240503750111. LE.DIFF)  GO  TO  9 
GO  TO  5 

G  IFCUU.LT.. 911312780288703)  GO  TO  8 

SECOND  NEARLY  LINEAR  DENSITY 

7  U=RANFC ISEED) 

H=RANFC ISEED) 

Z=U-U 

C  LET  U=  MAXCU.U)  AND  LET  W=MIWU.H) 

IFCU.GT.U)  GO  TO  101 
TEMP=U 
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u=w 

W=TEMP 

101  T=. 479727404222441+1. 105473BB1022070*W 
IF(U.LE. .87283497BB71790 )  GO  TO  9 

DIFF=EXP(-T*T*.5)/2.50S62827463100-.1800251910G85B3* 

*  (2.21B0358B716S471-ABSCT) ) 

IF( ABS(Z)*. 04926449S373128.LE. DIFF)  GO  TO  9 
GO  TO  7 

THIRD  NEARLY  LINEAR  DENSITY 

8  U=RANF( ISEED) 

W=RANF< ISEED) 

Z=U-W 

C  LET  U=  MAX(U*U)  AND  LET  H=MIN(U»U) 

IF(U.GT.U)  GO  TO  102 

TEMP=U 

U=U 

14= TEMP 

102  T=. 47972740422244 1- . 5955071380 1 S940*W 
IFCU.LE.. 805577924423817)  GO  TO  9 

DIFF=EXP( -T*T* . 5 ) ^2 . 50BB28274S3100-. 1800251910E85B3* 

*  (2.21S0358G716G471-ABSCT) ) 

IF  C  ABS ( Z ) * . 05337754950S88B . LE . DIFF )  GO  TO  9 
GO  TO  8 
9  X=T 

IFCZ.GE.0.0)  X=-X 

RETURN 

END 
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